Phase Transitions in Transportation Networks with Nonlinearities 
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We investigate a model of transportation networks with nonlinear elements which may represent 
local shortage of resources. Frustrations arise from competition for resources. When the initial 
resources are uniform, different regimes with discrete fractions of satisfied nodes are observed, re- 
sembling the Devil's staircase. We demonstrate how functional recursions are converted to simple 
recursions of probabilities. Behavior similar to those in the vertex cover or close packing problems 
are found. When the initial resources are bimodally distributed, increases in the fraction of rich 
nodes induce a glassy transition, entering an algorithmically hard regime. 

PACS numbers: 02.50.-r, 05.20.-y, 89.20.-a 



The study of currents and flows in networks is one 
of the most important problems in physics and many 
other areas of application Besides electric circuits 
transporting electric currents, there are transportation 
networks, communications networks, hydraulic networks, 
mammalian circulatory systems and vascular systems in 
plants @, OJ. A unified approach to these problems is 
facilitated by the minimization of cost functions. For ex- 
ample, they may represent the dissipation energy (via 
Thomson's principle for electric currents) [l| or time de- 
lays in communications networks. There is a close re- 
lation between the flow patterns and the cost functions. 
For example, it was found that the flow pattern is tree- 
like when the cost function is convex, but loopy other- 
wise [H] . These cost functions are continuous and often 
describe resistive flows. On the other hand, networks 
may contain nonlinear elements such as step-like penal- 
ties that may affect the flow distribution. 

The inclusion of nonlinear elements can drastically 
modify the flow patterns in the network. Nonlinearities 
are often represented by cost functions that are discon- 
tinuous in the currents. This creates numerous choices 
in deciding the current-carrying links and the idle ones. 
As shown in this paper, the flow patterns can demarcate 
clusters of nodes receiving currents. Clusters formed by 
similar energetic considerations have been found to play 
an important role in disordered systems such as the ran- 
dom field Ising model (RFIM) [4J, giving rise to the so- 
called Griffiths singularities and cascades of phase tran- 
sitions Q. As we shall see, similar behavior can also be 
found in nonlinear networks. 

The origin of these interesting phenomena can be 
traced the presence of frustrations, which refers to the 
conflict bewteen competing interaction energies in the 
system [6]. This connects such transportation networks 
with a broad class of network systems in which frustra- 
tion is inherent, including the /^-satisfiability Q, the 
graph coloring problem Q , and error-correcting codes @ . 
Spin glasses, the prototype of frustrated systems, provide 
the statistical mechanics to study these systems |10| . 

Transportation networks consist of nodes with either 
surplus or shortage of resources, and an important prob- 
lem is to distribute them to achieve a networkwide sat- 



isfaction with a minimum transportation cost [111 ]. This 
problem is important in load balancing in computer net- 
works and network flow of commodities [l2| . The focus 
of this paper is on the case where optimization of re- 
sources can proceed by penalizing nodes with shortages. 
In applications such as communications networks, such 
shortages can be detrimental to the performance. Hence 
it is interesting to consider step- like shortage costs, which 
give rise to unique behavior and physical picture absent 
in the previous models [TlT ] . Frustrations arise from com- 
petition for resources among connected nodes. Numer- 
ous metastable states emerge, leading to typical glassy 
behavior. 

Many NP-complete problems in computational com- 
plexity theory [l3j exhibit glassy behavior. As we shall 
see, the main results in this paper are that changes in 
the shortage per node induce phase transitions to glassy 
behavior. A new feature in our model is a cascade of 
phase transitions resembling the Devil's staircase in the 
glassy phase, and the first cascade is similar to the ver- 
tex cover problem [14j . As a crucial difference from other 
NP-complete problems, the present problem involve con- 
tinuous variables which complicate the analyses. We 
demonstrate how recursions of functions with continu- 
ous variables can be converted to simple recursions of 
probabilities. 

Specifically, we consider a dilute network of N nodes, 
labelled i — 1 . . . N. Each node is connected randomly 
to c neighbors with the symmetric connectivity matrix 
Aij = 1,0 for connected and unconnected node pairs re- 
spectively. Each node i has initial resource or capacity 
Aj, randomly drawn from a distribution of /d(A,). Posi- 
tive and negative values of Aj correspond to surplus and 
shortage of resources respectively. These resources can 
be transported through connected links. We denote by 
Uij = —yji the current of resources from node j to node 
i. The system minimizes the cost function 
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is the final resources of node i, and 
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Q(x) — 1 when x > 0, and otherwise. In the first term, 
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FIG. 2: The onset of different types of clusters of satisfied 
nodes for c = 3, with filled and unfilled circles representing 
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FIG. 1: (Color online) Simulation results of average energy 
per node, the fraction of unsatisfied nodes and the maximum 
cluster size of satisfied nodes after optimization, as a function 
of A/u. E' = 2E/u 2 . Parameters: c = 3, N = 100, 100 
samples and 1000 flips. Inset: the maximum cluster size of 
satisfied nodes shown with a larger vertical scale. 

each unsatisfied node yields a shortage cost of u 2 /2. The 
second term is the transportation cost of resources. 

We first look for phase transitions in the model by 
numerical simulations. To formulate an algorithm, we 
introduce a variable Si = ±1 for each node i and restrict 
its resources by s;£i > 0. Introducing Lagrange mul- 
tipliers ^ for the resource constraint, we minimize the 
Lagrangian 
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with the Kiihn- Tucker conditions /XiSj£i = and ^ < 
0. Optimizing L with respect to yij , one obtains yij = 
fijSj-fj,iSi and /x, = min[0, (A,-+£V A lj u j s j )/s l c]. Given 
a particular set of {s^}, we iterate these equations to find 
the corresponding set of {ui}. The set of optimal {si} is 
found by an approach similar to the the GSAT algorithm 
by comparing the Lagrangian in Eq. @ for each 
choice of {si}. 

We first consider networks with uniform capacity (A; = 
A for all i). As shown in Fig. [T]for c = 3, three phases 
can be identified: (1) unsatisfied phase for A/u < — v3, 
in which all the nodes are unsatisfied and E/N — u 2 /2; 

(2) partially satisfied phase for — yo < A/u < 0, in which 
only some nodes are satisfied, and < E/N < u 2 /2; 

(3) satisfied phase for A/u > 0, in which all nodes are 
satisfied and E/N = 0. 

Unlike the relatively smooth decreasing trend of en- 
ergy, the fraction of unsatisfied nodes is a discontinuous 
function of A/u, showing abrupt jumps at threshold val- 
ues of A/u. The step size of the curve decreases as A/u 
increases, and gradually becomes unresolvable by the nu- 
merical experiments. This resembles the Devil's staircase 
observed in the circle map and other dynamical systems 
[l6l |. These threshold values of A/u mark the position 



FIG. 3: (Color online) (a) A typical set of cavity energy 
functions E k (y) at c = 3, A/u = -1.22 > -y/J/2. (b- 
d) The closed set of Ek(y) at c = 3, A/u — —5/3 with 
7 = u 2 /2 — A/2c, corresponding to (b) the u-state, (c) the 
s-state and (d) the b-state. 



at which certain configurations of the satisfied nodes be- 
come energetically stable. 

To further confirm this picture, we measure the aver- 
age maximum cluster size of satisfied nodes in the sam- 
ples. Abrupt jumps of the cluster size are observed at the 
same threshold values as shown in Fig. [1] This indicates 
that new types of clusters of satisfied nodes are formed at 
each jump, as shown in Fig. [2]for c = 3. The situation is 
similar to the cascades of phase transitions in RFIM due 
to the onsets of different clusters ||. As shown in the 
inset of Fig. [2 the maximum cluster size increases dras- 
tically to O(N) around A/u « —0.5. This corresponds 
to a percolation-like transition at which isolated clusters 
of satisfied nodes become connected. 

We carry out the analysis using the Bethe approxima- 
tion, since the networks have a tree-like structure locally. 
Messages are passed from the vertices of the sub-trees 
to the next layer, resulting in a recursion relation of the 
messages. In the cavity approach, these messages are the 
cavity energy functions, denoted by Ej{yij) for the en- 
ergy of the tree terminated at a link from vertex j to its 
ancestor i, when a current yij is drawn from j to i [Til ] . 
They can be expressed in terms of the energies of its de- 
scendents k = 1, . . . , c— 1, 
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where £, j {y lJ ) = Aj + £ fc y jk - Vi] . 

In general, the iteration of Eq. ([3]) results in a dis- 
tribution of the cavity energy functions, as shown in 
Fig. Ela). However, in the regime —\[c < A/u < 
— yjc(c— l)/(c+ 1), it can be shown analytically that 
there are c functional forms of E k (y) forming a closed set 
of solutions to the recursion relation in Eq. The three 
functions for c = 3, referred to as the u-, s- and 6-states, 



are shown in Fig. (3^b-d). Respectively, they correspond 
to states favoring unsatisfaction, satisfaction, and testa- 
bility, and have absolute minima at y = 0, y — A/c, and 
both y = and A/c. Furthermore, numerical iterations 
of Eq. ([3]) starting with random Ek (y) show that this set 
of solutions is stable. For c > 3, the closed set consists 
of more states having absolute minima at y = 0, and lo- 
cal minima at different energies at y = A/c, similar to 
the it-state in Fig. [3l As we shall see, the closure prop- 
erty of the c states greatly simplifies Eq. With the 
u- and 6-states denoted as the U -state, and the s-state 
by the S-state, their recursion relations are given by an 
'^[/-recursion" : assign a vertex to an S-state if all its 
c — 1 descendents are in the [/-state, and to a [/-state 
otherwise. 

The recursion rules imply that optimization is achieved 
by the close packing of satisfied nodes in the network, 
with the constraint that they do not form clusters. Hence 
we call the regime — -Jc < A/u < — yc(c — l)/(c + 1) 
the single-sat regime. This reduces the problem to the 
vertex cover problem [bj ] . Since there is at most one sat- 
isfied node at the end of each link, the maximization of 
the number of satisfied nodes is equivalent to the mini- 
mization of the covered set size in the vertex cover prob- 
lem. Alternatively, the model can be considered as the 
close packing limit of a lattice glass model 17] . Both 
models exhibit glassy behavior, and phases with multi- 
ple metastable states are found therein, corresponding to 
the computationally hard phases. 

The description of the glassy behavior can be ap- 
proached by the replica analysis In the so-called 
replica symmetric (RS) ansatz, one assumes that the net- 
work behavior is dominated by a single ground state. 
However, we find that this ansatz is not stable in the 
entire single-sat phase. Instead, the network behavior is 
described by the so-called full replica symmetry-breaking 
(FRSB) ansatz, which assumes an infinite hierarchy of 
metastable states. Indeed, we have computed the opti- 
mized energy in the one-step RSB (1RSB) approxima- 
tion, assuming that the optimized states with energy E 
are distributed exponentially according to exp[NY,(E)], 
where £(£/) is the so-called configurational entropy [la |. 
Futhermore, the 1RSB ansatz predicts that 'S(E) exists 
up to an energy Ed above the ground state energy E s , and 
the numerous metastable states prevent practical algo- 
rithms from converging to the true ground state, result- 
ing in a dynamical transition. Specifically, we find that 
in the single-sat regime, E/N = A 2 /6 + f u {u 2 /2 - A 2 /6), 
where f u is the fraction of unsatisfied nodes. In the RS 
approximation, /„ = 0.545, whereas in the 1RSB approx- 
imation, f u — 0.549 and 0.550 at its E s and Ed respec- 
tively. Note that f u at the dynamical transition is close 
to the simulation result of f u — 0.551. 

When A increases above the single-sat regime, clusters 
of two satisfied nodes appear. This double-sat regime ex- 
ists in the range —a/3/2 < A/u < — -»/21/25 for c = 3. 
This is similar to the close packing limit of the Bethe glass 
model, in which each occupied site can have at most one 
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FIG. 4: (a) A five-node configuration composing of 3 isolated 
satisfied nodes (filled circles). It is stable when A/u < — 1. (b) 
A five-node configuration composing of 2 two-node clusters. 
It replaces the configuration in (a) when A/u > — 1. 

occupied neighbor. However, the present model is richer 
in behavior, as indicated by the jumps in the fraction 
of unsatisfied nodes in the double-sat regime of Fig. [T] 
They mark the positions at which configurations increas- 
ingly dominated by two-node clusters become energeti- 
cally stable, when A increases. An example of such an 
energy switch is shown in Fig. 2] Indeed, such thresh- 
old values are expected at A/u = — \/ (3ra — 3)/ (2n — 1) 
where n = 2,3, ... . Only two of those thresholds are 
visible in Fig. [TJ probably due to the absence of config- 
urations for larger n in networks with N = 100, and the 
limitations of the search algorithm. Note that the double- 
sat regime may exhibit behaviors described by different 
RSB ansatz as A/u varies. It is also interesting to study 
the possible change of RSB behavior as A/u further in- 
creases up the Devil's staircase. 

Next, we consider a simple case with quenched disorder 
of the node capacities. In this case, Aj is drawn from a 
bimodal distribution, namely, Aj = A± with probabilities 
f±, where A_/it < -\/3 < A + /u < -y/3/2 for c = 
3, and / + + /_ = 1. A± are referred to as the rich 
and poor nodes respectively. The recursion relations for 
the rich nodes follow the "^[/-recursion" , whereas the 
poor nodes are always in the [/-state. Note that the end 
points of this range are /+ = and 1, corresponding to 
the unsatisfied and partially satisfied phases respectively. 
Since these two phases are in the RS and RSB regimes 
respectively, we expect that there is a phase transition 
when / + increases from to 1. 

Applying the RS ansatz to the region of low / + , we let 
Pu,k be the probability that node k is in the [/-state. Its 
recursion relation can be written as 



Pu,j = ^A 3 -,A_ + <5a 3 , 
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which corresponds to an intense simplificaton of Eq. ([3]) 
In the easy regime at low / + , a stable fixed point solution 
with all p Ut k = or 1 exists. The site average (p u ) can 
thus be obtained from the equation (p u ) = /_ + — 
(Pm) c_1 )- In terms of algorithms, optimal network states 
in this regime can be obtained by the so-called belief 
propagation (BP) algorithm, initializing the messages to 
or 1. 

The stability of the RS solution can be studied by con- 
sidering the propagation of fluctuations {{5p u .k) 2 ) under 
the recursion relation Eq. ([J) fl9| . This leads to the 
Almeida-Thouless (AT) stability condition, which reads 
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(a) The fraction of nodes with p Ut k = 0, 1, and < 
1 for c = 3. (b) The phase diagram. Inset: Simulation 
of the fraction of messages from rich nodes with BP 
ence on dilute networks with c = 3. 



In the RS regime, (pi) = (p u ) since p Ut k = or 1. This 



implies an AT transition at /+ T 
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As shown in Fig. EJa), free nodes with < p u ,k < 1 
start to exist when / + > f+ T , analogous to the vertex 
cover [3] and graph coloring problems Q. This charac- 
terizes the hard region with multiple states. 

As shown in the /+-c space in Fig. [5Jb), the easy and 
hard phases exist below and above the AT line respec- 
tively. In the large c limit, approaches e/c. This 
result has an interesting connection with the vertex cover 
problem. Considering the cover set as the set of unsat- 
isfied nodes, all links involving poor nodes are covered. 
The remaining links are those among the rich nodes, and 
the problem of minimizing the cover set size reduces to 
one that minimizes the subset size of covered nodes in the 



subnetwork of rich nodes. In the large c limit, this sub- 
network has a Poissonian connectivity distribution with 
an average c/ + . The result c/^ T = e agrees with the 
point of RS instability derived in [l4| • 

The AT transition has implications on the algorithms. 
We consider the BP algorithm initialized with messages 
and 1. As shown in Fig. [5fb) inset, effectively all mes- 
sages from rich nodes converge in the easy regime. How- 
ever, a significant fraction of messages fluctuates between 
and 1 above /+ T , indicating the breakdown of the RS 
ansatz. Algorithmically, decimation procedures, such as 
those used in the survey propagation (SP) algorithm 0], 
are required. 

In summary, we have studied how nonlinearities af- 
fect the flow patterns in transportation networks. In 
the case of uniform capacity, phase transitions resembling 
the Devil's staircase reveal the cascades of clustered flow 
patterns. In the single-sat regime with a closed set of 
only a few cavity energy functions, the flow pattern has 
a correspondence with the vertex cover problem. Glassi- 
ness arises from the frustration in competitions for re- 
sources. In the case with quenched disorder of the ca- 
pacities, an increase in the fraction of rich nodes induces 
a phase transition from an easy phase to a hard one, 
in which message-passing algorithms experience conver- 
gence problems. These features are relevant to general 
network optimization problems with nonlinear elements. 
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